1 Setup

1.1 Load functions and set directories

Set knitting options

# global knitting options for automatic saving of all plots as .png and .pdf. Also sets cache directory.
knitr::opts_chunk$set(
  dev = c("png", "pdf"), fig.keep = "all",
  dev.args = list(pdf = list(encoding = "WinAnsi", useDingbats = FALSE)),
  fig.path = file.path("fig_output/", paste0(gsub("\\.[Rr]md", "/", knitr::current_input()))),
  cache.path = file.path("cache/", paste0(gsub("\\.[Rr]md", "/", knitr::current_input())))
)

Load packages

library(dada2); packageVersion("dada2") # sequence filtering and sample inference with dada2
## [1] '1.16.0'
library(ShortRead); packageVersion("ShortRead") # package for loading, displaying, and handling sequences
## [1] '1.46.0'
library(tidyverse); packageVersion("tidyverse") # handy piping and data operations
## [1] '1.3.0'
library(rjson); packageVersion("rjson") # read json output from figaro
## [1] '0.2.20'
# source all relevant scripting files
# paths.R contains paths to local data and programs
source(file.path("scripts", "paths.R"))

Set directories

# Set up names of sub directories to stay organized
preprocess_path <- file.path(project_path_OM18, "01_preprocess")
    demultiplex_path <- file.path(preprocess_path, "demultiplexed")
        demultiplex_stat_summ_path <- file.path(demultiplex_path, "stat_summ") # statistical summary output from idemp
        demultiplex_unassigned_path <- file.path(demultiplex_path, "unassigned")
        demultiplex_assigned_path <- file.path(demultiplex_path, "assigned")
    figaro_output_path <- file.path(preprocess_path, "Figaro_output")
filter_path <- file.path(project_path_OM18, "02_filter")

1.2 Organize and rename demultiplexing output files

(demux_output_raw_names <- sort(list.files(demultiplex_path, full.names = FALSE)))
## [1] "assigned"

Move demultiplexing output files into different directories depending on whether they are statistical summaries, unassigned fastq files, or assigned fastq files

# Create directory to hold idemp statistical summary files
if (!dir.exists(demultiplex_stat_summ_path)) dir.create(demultiplex_stat_summ_path)
# Now, rename (move) files
file.rename(from = list.files(demultiplex_path, pattern="decode", full.names = TRUE), to = file.path(demultiplex_stat_summ_path, list.files(demultiplex_path , pattern="decode", full.names = FALSE)))
## logical(0)
# Create directory to hold idemp unassigned fastq files
if (!dir.exists(demultiplex_unassigned_path)) dir.create(demultiplex_unassigned_path)
# Now, rename (move) files
file.rename(from = list.files(demultiplex_path, pattern="unsigned", full.names = TRUE), to = file.path(demultiplex_unassigned_path, list.files(demultiplex_path , pattern="unsigned", full.names = FALSE)))
## logical(0)
# Create directory to hold idemp assigned fastq files
if (!dir.exists(demultiplex_assigned_path)) dir.create(demultiplex_assigned_path)
# Now, rename (move) files
file.rename(from = list.files(demultiplex_path, pattern="fastq", full.names = TRUE), to = file.path(demultiplex_assigned_path, list.files(demultiplex_path , pattern="fastq", full.names = FALSE)))
## logical(0)

Now we read in the names of the demultiplexed and assigned fastq files.

fnFs_raw_names <- sort(list.files(demultiplex_assigned_path, pattern = "R1", full.names = TRUE))
fnRs_raw_names <- sort(list.files(demultiplex_assigned_path, pattern = "R2", full.names = TRUE))

Extract sample names for later use

# Extract sample names
sample_names_F <- fnFs_raw_names %>% basename() %>% str_remove("Undetermined_S0_L001_R1_001.fastq.gz_") %>% str_remove(".fastq.gz")

sample_names_R <- fnRs_raw_names %>% basename() %>% str_remove("Undetermined_S0_L001_R2_001.fastq.gz_") %>% str_remove(".fastq.gz")

# Check that forward and reverse file names match
if(!identical(sample_names_F, sample_names_R)) stop("Forward and reverse file names do not match.")

# print first few names
sample_names_F %>% head()
## [1] "BA1A100_10_2J"  "BA1A100_22_2D"  "BA1A100_45_2C"  "BA1A5566_10_1D"
## [5] "BA1A5566_22_1C" "BA1A5566_45_1J"
# Repeat the above with reverse reads
R1_names_incremented <- vector("character", length(sample_names_F))
for (i in seq_along(R1_names_incremented )) {
  R1_names_incremented[[i]] <- paste0(sample_names_F[[i]], "_S", i, "_L001_R1_001.fastq.gz")
}

# Now, rename files
file.rename(from = fnFs_raw_names, to = file.path(demultiplex_assigned_path, R1_names_incremented))
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
fnFs <- sort(list.files(demultiplex_assigned_path, pattern = "R1", full.names = TRUE)) # re-assign full names of forward reads, now with Illumina naming convention for use with Figaro.

# Print first few lines of fnFs
fnFs %>% head()
## [1] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/demultiplexed/assigned/BA1A100_10_2J_S1_L001_R1_001.fastq.gz" 
## [2] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/demultiplexed/assigned/BA1A100_22_2D_S2_L001_R1_001.fastq.gz" 
## [3] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/demultiplexed/assigned/BA1A100_45_2C_S3_L001_R1_001.fastq.gz" 
## [4] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/demultiplexed/assigned/BA1A5566_10_1D_S4_L001_R1_001.fastq.gz"
## [5] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/demultiplexed/assigned/BA1A5566_22_1C_S5_L001_R1_001.fastq.gz"
## [6] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/demultiplexed/assigned/BA1A5566_45_1J_S6_L001_R1_001.fastq.gz"
# Repeat the above with reverse reads
R2_names_incremented <- vector("character", length(sample_names_R))
for (i in seq_along(R2_names_incremented )) {
  R2_names_incremented[[i]] <- paste0(sample_names_R[[i]], "_S", i, "_L001_R2_001.fastq.gz")
}

file.rename(from = fnRs_raw_names, to = file.path(demultiplex_assigned_path, R2_names_incremented))
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
fnRs <- sort(list.files(demultiplex_assigned_path, pattern = "R2", full.names = TRUE))

Count and display reads for input

# count forward reads
input_read_count_F <- countLines(demultiplex_assigned_path, pattern="_R1_001.fastq")/4

# make forward read counts into tbl with sample name
input_read_count_F_tbl <- input_read_count_F %>% as_tibble(rownames = "file_name") %>% rename(reads_F = value) %>% mutate(sample_id = file_name %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R1_001.fastq.gz"))

# repeat for reverse

input_read_count_R <- countLines(demultiplex_assigned_path, pattern="_R2_001.fastq")/4

input_read_count_R_tbl <- input_read_count_R %>% as_tibble(rownames = "file_name") %>% rename(reads_R = value) %>% mutate(sample_id = file_name %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R2_001.fastq.gz"))


# Join forward and reverse together into one tbl
input_read_count_FR <- input_read_count_F_tbl %>% select(-file_name) %>% left_join(input_read_count_R_tbl %>% select(-file_name), by = "sample_id") %>% select(sample_id, everything())

# print. reads counts are same in forward and reverse, which is what we want.
input_read_count_FR

Plot read counts

plot_read_count_input <- input_read_count_FR %>% ggplot(aes(
  x = fct_reorder(sample_id, desc(reads_F)),
  y = reads_F
)) +
  geom_bar(stat = "identity") +
  scale_y_continuous(name = "Forward reads") +
  scale_x_discrete(name = "Sample ID") +
  theme_classic(base_size = 9)+
  theme(
    axis.text.x = element_text(angle = 90, vjust = .5, hjust = 1),
    legend.position = "bottom"
  )

plot_read_count_input

2 Check primer sequences

The 515 (forward) (Parada) and 806R (reverse) (Apprill) primers (EarthMicrobiome) were used to amplify this dataset. We record the DNA sequences for those primers.

However, the sequencing strategy used here is not supposed to sequence the primers.

FWD <- "GTGYCAGCMGCCGCGGTAA" # 515F (Parada)
REV <- "GGACTACNVGGGTWTCTAAT" # 806R (Apprill)

(FWD_length <- str_length(FWD))
## [1] 19
(REV_length <- str_length(REV))
## [1] 20

That said, let’s check for any primer sequence matches in the sequences.

allOrients <- function(primer) {
    # Create all orientations of the input sequence
    require(Biostrings)
    dna <- DNAString(primer)  # The Biostrings works w/ DNAString objects rather than character vectors
    orients <- c(Forward = dna, Complement = complement(dna), Reverse = reverse(dna), 
        RevComp = reverseComplement(dna))
    return(sapply(orients, toString))  # Convert back to character vector
}
FWD.orients <- allOrients(FWD)
REV.orients <- allOrients(REV)
FWD.orients
##               Forward            Complement               Reverse 
## "GTGYCAGCMGCCGCGGTAA" "CACRGTCGKCGGCGCCATT" "AATGGCGCCGMCGACYGTG" 
##               RevComp 
## "TTACCGCGGCKGCTGRCAC"

We are now ready to count the number of times the primers appear in the forward and reverse read, while considering all possible primer orientations. Identifying and counting the primers on one set of paired end FASTQ files is sufficient, assuming all the files were created using the same library preparation, so we’ll just process the first sample.

# write function to check for primer sequences
primerHits <- function(primer, fn) { # Counts number of reads in which the primer is found
    nhits <- vcountPattern(primer, sread(readFastq(fn)), fixed = FALSE)
    return(sum(nhits > 0))
}
# print file name
print(paste("file_name: ", basename(fnFs[[1]])))
## [1] "file_name:  BA1A100_10_2J_S1_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs[[1]])/4))
## [1] "total reads: 28746"
# check for primers
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs[[1]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnRs[[1]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs[[1]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnRs[[1]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       7
## REV.ForwardReads       0          0       0       9
## REV.ReverseReads       0          0       0       0
# print file name
print(paste("file_name: ", basename(fnFs[[2]])))
## [1] "file_name:  BA1A100_22_2D_S2_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs[[2]])/4))
## [1] "total reads: 40899"
# check for primers
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs[[2]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnRs[[2]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs[[2]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnRs[[2]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       1
## REV.ForwardReads       0          0       0       2
## REV.ReverseReads       0          0       0       0
# print file name
print(paste("file_name: ", basename(fnFs[[3]])))
## [1] "file_name:  BA1A100_45_2C_S3_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs[[3]])/4))
## [1] "total reads: 31399"
# check for primers
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs[[3]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnRs[[3]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs[[3]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnRs[[3]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       1
## REV.ForwardReads       0          0       0       3
## REV.ReverseReads       1          0       0       0

There were a few (<10) detections of FWD primer sequences in the reverse reads in its reverse-complement orientation, and vice verse for the REV primer. We will check if primer sequences are still found later after we truncate lower quality reads. So now let’s inspect read quality profiles.

3 Inspect read quality profiles

It’s important to get a feel for the quality of the data that we are using. To do this, we will plot the quality of some of the samples.

From the dada2 tutorial: >In gray-scale is a heat map of the frequency of each quality score at each base position. The median quality score at each position is shown by the green line, and the quartiles of the quality score distribution by the orange lines. The red line shows the scaled proportion of reads that extend to at least that position (this is more useful for other sequencing technologies, as Illumina reads are typically all the same lenghth, hence the flat red line).

# Forward reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnFs) <= 20) {
  plotQualityProfile(fnFs)
} else {
  rand_samples <- sample(size = 20, 1:length(fnFs)) # grab 20 random samples to plot
  plotQualityProfile(fnFs[rand_samples])
}

# Reverse reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnRs) <= 20) {
  plotQualityProfile(fnRs)
} else {
  rand_samples <- sample(size = 20, 1:length(fnRs)) # grab 20 random samples to plot
  plotQualityProfile(fnRs[rand_samples])
}

4 Filter parameter optimization with Figaro

Create a directory for Figaro output

# Create directory to hold the output from figaro
if (!dir.exists(figaro_output_path)) dir.create(figaro_output_path)

I couldn’t get figaro (v1.1.1) to run by system2() in R due to a shell permission denied error, but it worked on the command line. This chunk generates the exact text I entered in the command line after changing to the figaro directory. I then just read in the figaro output JSON and continue from there.

# text to change directory in terminal to figaro 
paste("cd", figaro_path) %>% str_remove("/figaro.py") %>% cat()
## cd /Users/melo.d/opt/figaro
# primers were not sequenced
# sequences are 1X150 bp
# Tried 240, lost nearly all reads at merge step.
# Try again with 260
figaro_amplicon_length <- 260

figaro_min_overlap <- 20

figaro_output_file_name <- "Figaro_params_OM18.json"

# using 1 for primer lengths because even though primers were not sequenced, the current version of figaro does not allow 0 as a primer length.

flags_figaro <- paste("--ampliconLength", figaro_amplicon_length, "--forwardPrimerLength", 1, "--reversePrimerLength", 1, "--outputFileName", figaro_output_file_name, "--inputDirectory", demultiplex_assigned_path, "--outputDirectory", figaro_output_path, "--minimumOverlap", figaro_min_overlap)

paste("python3", "figaro.py", flags_figaro) %>% cat()
## python3 figaro.py --ampliconLength 260 --forwardPrimerLength 1 --reversePrimerLength 1 --outputFileName Figaro_params_OM18.json --inputDirectory /Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/demultiplexed/assigned --outputDirectory /Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/01_preprocess/Figaro_output --minimumOverlap 20

Read in Figaro JSON output, top 5 suggested trim parameters.

# Give the input file name to the function.
figaro_result <- fromJSON(file = file.path(figaro_output_path, figaro_output_file_name))

figaro_result[1:5] %>% paste(collapse = '\n') %>% cat()
## list(trimPosition = c(147, 135), maxExpectedError = c(2, 2), readRetentionPercent = 91.72, score = 89.7169714831002)
## list(trimPosition = c(148, 134), maxExpectedError = c(2, 2), readRetentionPercent = 91.71, score = 89.7116993626006)
## list(trimPosition = c(146, 136), maxExpectedError = c(2, 2), readRetentionPercent = 91.7, score = 89.7022095457014)
## list(trimPosition = c(149, 133), maxExpectedError = c(2, 2), readRetentionPercent = 91.7, score = 89.7022095457014)
## list(trimPosition = c(145, 137), maxExpectedError = c(2, 2), readRetentionPercent = 91.69, score = 89.6853387601027)
# Save and print Figaro-suggested parameters for dada2::filterAndTrim
(truncLen_best_figaro <- figaro_result[[1]]$trimPosition)
## [1] 147 135
(maxEE_best_figaro <- figaro_result[[1]]$maxExpectedError)
## [1] 2 2

5 Filter and trim for quality

5.1 Organize file paths and directories

# set file path for quality-filtered read
filter_path <- file.path(project_path_OM18, "02_filter") 

# Put filtered reads into separate sub-directories for big data workflow

if (!dir.exists(filter_path)) dir.create(filter_path)
    subF_path <- file.path(filter_path, "preprocessed_F") 
    subR_path <- file.path(filter_path, "preprocessed_R") 
if (!dir.exists(subF_path)) dir.create(subF_path)
if (!dir.exists(subR_path)) dir.create(subR_path)
    
# Move demultiplexed and assigned R1 and R2 to separate forward/reverse sub-directories

fnFs.Q <- file.path(subF_path,  basename(list.files(demultiplex_assigned_path , pattern="R1", full.names = TRUE))) 
fnRs.Q <- file.path(subR_path,  basename(list.files(demultiplex_assigned_path , pattern="R2", full.names = TRUE)))

file.rename(from = list.files(demultiplex_assigned_path , pattern="R1", full.names = TRUE), to = fnFs.Q)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
file.rename(from = list.files(demultiplex_assigned_path , pattern="R2", full.names = TRUE), to = fnRs.Q)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# File parsing; create file names and make sure that forward and reverse files match
filtpathF <- file.path(subF_path, "filtered") # files go into preprocessed_F/filtered/
filtpathR <- file.path(subR_path, "filtered") # ...
fastqFs <- sort(list.files(subF_path, pattern="fastq.gz"))
fastqRs <- sort(list.files(subR_path, pattern="fastq.gz"))
if(length(fastqFs) != length(fastqRs)) stop("Forward and reverse files do not match.")

Before chosing sequence variants, we want to trim reads where their quality scores begin to drop (the truncLen and truncQ values) and remove any low-quality reads that are left over after we have finished trimming (the maxEE value).

5.2 Filter and trim

# No argument pasted to trimLeft because the primers were not sequenced
filt_out <- filterAndTrim(fwd = file.path(subF_path, fastqFs), filt = file.path(filtpathF, fastqFs), rev = file.path(subR_path, fastqRs), filt.rev = file.path(filtpathR, fastqRs), truncLen = truncLen_best_figaro, maxEE = maxEE_best_figaro, truncQ = 2, maxN = 0, rm.phix = TRUE, compress = TRUE, verbose = TRUE, multithread = TRUE)

# look at how many reads were kept
head(filt_out)
##                                        reads.in reads.out
## BA1A100_10_2J_S1_L001_R1_001.fastq.gz     28746     27391
## BA1A100_22_2D_S2_L001_R1_001.fastq.gz     40899     39238
## BA1A100_45_2C_S3_L001_R1_001.fastq.gz     31399     30259
## BA1A5566_10_1D_S4_L001_R1_001.fastq.gz    51141     48913
## BA1A5566_22_1C_S5_L001_R1_001.fastq.gz    42788     40767
## BA1A5566_45_1J_S6_L001_R1_001.fastq.gz    31469     29943
# summary of samples in filt_out by percentage
filt_out %>% 
  data.frame() %>% 
  mutate(Samples = rownames(.),
         percent_kept = 100*(reads.out/reads.in)) %>%
  select(Samples, everything()) %>%
  summarise(min_remaining = paste0(round(min(percent_kept), 2), "%"), 
            median_remaining = paste0(round(median(percent_kept), 2), "%"),
            mean_remaining = paste0(round(mean(percent_kept), 2), "%"), 
            max_remaining = paste0(round(max(percent_kept), 2), "%"))

5.3 Check if primers were removed by filterAndTrim

Generate sorted file name vector for filtered files

fnFs_filt <- sort(list.files(filtpathF, full.names = TRUE))

# print first 6 lines
fnFs_filt %>% head()
## [1] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/02_filter/preprocessed_F/filtered/BA1A100_10_2J_S1_L001_R1_001.fastq.gz" 
## [2] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/02_filter/preprocessed_F/filtered/BA1A100_22_2D_S2_L001_R1_001.fastq.gz" 
## [3] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/02_filter/preprocessed_F/filtered/BA1A100_45_2C_S3_L001_R1_001.fastq.gz" 
## [4] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/02_filter/preprocessed_F/filtered/BA1A5566_10_1D_S4_L001_R1_001.fastq.gz"
## [5] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/02_filter/preprocessed_F/filtered/BA1A5566_22_1C_S5_L001_R1_001.fastq.gz"
## [6] "/Users/melo.d/Desktop/Research/Oman_fieldwork/Samail_16S_compilation_fq_outputs/2018/02_filter/preprocessed_F/filtered/BA1A5566_45_1J_S6_L001_R1_001.fastq.gz"
fnRs_filt <- sort(list.files(filtpathR, full.names = TRUE))
# print file name
print(paste("file_name: ", basename(fnFs_filt[[1]])))
## [1] "file_name:  BA1A100_10_2J_S1_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs_filt[[1]])/4))
## [1] "total reads: 27391"
# tally primer hits in file
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[1]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[1]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[1]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[1]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       0
## REV.ForwardReads       0          0       0       2
## REV.ReverseReads       0          0       0       2
# print file name
print(paste("file_name: ", basename(fnFs_filt[[2]])))
## [1] "file_name:  BA1A100_22_2D_S2_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs_filt[[2]])/4))
## [1] "total reads: 39238"
# tally primer hits in file
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[2]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[2]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[2]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[2]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       0
## REV.ForwardReads       0          0       0       0
## REV.ReverseReads       0          0       0       0
# print file name
print(paste("file_name: ", basename(fnFs_filt[[3]])))
## [1] "file_name:  BA1A100_45_2C_S3_L001_R1_001.fastq.gz"
print(paste("total reads:", countLines(fnFs_filt[[3]])/4))
## [1] "total reads: 30259"
# tally primer hits in file
rbind(FWD.ForwardReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[3]]),
      FWD.ReverseReads = sapply(FWD.orients, primerHits, fn=fnFs_filt[[3]]),
      REV.ForwardReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[3]]),
      REV.ReverseReads = sapply(REV.orients, primerHits, fn=fnFs_filt[[3]]))
##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       0
## REV.ForwardReads       0          0       0       1
## REV.ReverseReads       0          0       0       1

Primer sequences appear in ~1 in every 10 000 reads. At such low primer counts, additional filtering of primers is probably not needed, and won’t be performed here.

5.4 Plot the quality of the filtered and cut fastq files.

# Forward reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnFs_filt) <= 20) {
  plotQualityProfile(fnFs_filt)
} else {
  rand_samples <- sample(size = 20, 1:length(fnFs_filt)) # grab 20 random samples to plot
  plotQualityProfile(fnFs_filt[rand_samples])
}

# Forward reads
# If the number of samples is 20 or less, plot them all, otherwise, just plot 20 randomly selected samples
if( length(fnFs_filt) <= 20) {
  plotQualityProfile(fnFs_filt)
} else {
  rand_samples <- sample(size = 20, 1:length(fnFs_filt)) # grab 20 random samples to plot
  plotQualityProfile(fnFs_filt[rand_samples])
}

6 Infer sequence variants

In this part of the pipeline, dada2 will learn to distinguish error from biological differences using a subset of our data as a training set. After it understands the error rates, we will reduce the size of the dataset by combining all identical sequence reads into “unique sequences”. Then, using the dereplicated data and error rates, dada2 will infer the sequence variants in our data. Finally, we will merge the corresponding forward and reverse reads to create a list of the fully denoised sequences and create a sequence table from the result.

6.1 Learn the error rates

set.seed(100) # set seed to ensure that randomized steps are reproducible

# Learn forward error rates
errF <- learnErrors(fnFs_filt, nbases=1e8, multithread=TRUE)
## 106221906 total bases in 722598 reads from 44 samples will be used for learning the error rates.
# Learn reverse error rates
errR <- learnErrors(fnRs_filt, nbases=1e8, multithread=TRUE)
## 102431385 total bases in 758751 reads from 45 samples will be used for learning the error rates.

6.1.1 Plot Error Rates

The error rates for each possible transition (A→C, A→G, …) are shown. Points are the observed error rates for each consensus quality score. The black line shows the estimated error rates after convergence of the machine-learning algorithm. The red line shows the error rates expected under the nominal definition of the Q-score. Here, we want to see estimated error rates (black line) that are a good fit to the observed rates (points), and that the error rates drop with increased quality as expected. If these criteria aren’t observed, then it may be a good idea to increase the nbases parameter. This allows the machine learning algorithm to train on a larger portion of your data and may help improve the fit.

plotErrors(errF, nominalQ=TRUE)
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

plotErrors(errR, nominalQ=TRUE)
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

6.2 Dereplication, sequence inference, and merging of paired-end reads

Extract sample names and add them to vector of filtered file names

# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sample_names_F_filt <- fnFs_filt %>% basename() %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R1_001.fastq.gz")
sample_names_R_filt <- fnRs_filt %>% basename() %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R2_001.fastq.gz")

# Double check
if(!identical(sample_names_F_filt, sample_names_R_filt)) stop("Forward and reverse files do not match.")

names(fnFs_filt) <- sample_names_F_filt
names(fnRs_filt) <- sample_names_R_filt
# make lists to hold the loop output
mergers <- vector("list", length(sample_names_F))
names(mergers) <- sample_names_F
ddF <- vector("list", length(sample_names_F))
names(ddF) <- sample_names_F
ddR <- vector("list", length(sample_names_F))
names(ddR) <- sample_names_F

# For each sample, get a list of merged and denoised sequences
for(sam in sample_names_F) {
    cat("Processing:", sam, "\n")
    # Dereplicate forward reads
    derepF <- derepFastq(fnFs_filt[[sam]])
    # Infer sequences for forward reads
    ddF[[sam]] <- dada(derepF, err=errF, multithread=TRUE)
    # Dereplicate reverse reads
    derepR <- derepFastq(fnRs_filt[[sam]])
    # Infer sequences for reverse reads
    ddR[[sam]] <- dada(derepR, err=errR, multithread=TRUE)
    # Merge reads together
    merger <- mergePairs(ddF[[sam]], derepF, ddR[[sam]], derepR)
    mergers[[sam]] <- merger
}
## Processing: BA1A100_10_2J 
## Sample 1 - 27391 reads in 5354 unique sequences.
## Sample 1 - 27391 reads in 4332 unique sequences.
## Processing: BA1A100_22_2D 
## Sample 1 - 39238 reads in 9156 unique sequences.
## Sample 1 - 39238 reads in 7628 unique sequences.
## Processing: BA1A100_45_2C 
## Sample 1 - 30259 reads in 6639 unique sequences.
## Sample 1 - 30259 reads in 5308 unique sequences.
## Processing: BA1A5566_10_1D 
## Sample 1 - 48913 reads in 12046 unique sequences.
## Sample 1 - 48913 reads in 10884 unique sequences.
## Processing: BA1A5566_22_1C 
## Sample 1 - 40767 reads in 9901 unique sequences.
## Sample 1 - 40767 reads in 8388 unique sequences.
## Processing: BA1A5566_45_1J 
## Sample 1 - 29943 reads in 7540 unique sequences.
## Sample 1 - 29943 reads in 6008 unique sequences.
## Processing: CM2A_10_9C 
## Sample 1 - 33380 reads in 5257 unique sequences.
## Sample 1 - 33380 reads in 4124 unique sequences.
## Processing: CM2A_22_9J 
## Sample 1 - 29956 reads in 5612 unique sequences.
## Sample 1 - 29956 reads in 4498 unique sequences.
## Processing: CM2A_45_9D 
## Sample 1 - 35762 reads in 6843 unique sequences.
## Sample 1 - 35762 reads in 5854 unique sequences.
## Processing: NSHQ14_10_7D 
## Sample 1 - 24302 reads in 5745 unique sequences.
## Sample 1 - 24302 reads in 4870 unique sequences.
## Processing: NSHQ14_22_7C 
## Sample 1 - 23989 reads in 4909 unique sequences.
## Sample 1 - 23989 reads in 4174 unique sequences.
## Processing: NSHQ14_45_7J 
## Sample 1 - 25565 reads in 5570 unique sequences.
## Sample 1 - 25565 reads in 4903 unique sequences.
## Processing: NSHQ4_S_11D 
## Sample 1 - 34607 reads in 9635 unique sequences.
## Sample 1 - 34607 reads in 6683 unique sequences.
## Processing: NTC_1 
## Sample 1 - 1150 reads in 776 unique sequences.
## Sample 1 - 1150 reads in 732 unique sequences.
## Processing: NTC_10 
## Sample 1 - 2284 reads in 1386 unique sequences.
## Sample 1 - 2284 reads in 1246 unique sequences.
## Processing: NTC_11 
## Sample 1 - 2284 reads in 542 unique sequences.
## Sample 1 - 2284 reads in 476 unique sequences.
## Processing: NTC_12 
## Sample 1 - 2286 reads in 1092 unique sequences.
## Sample 1 - 2286 reads in 953 unique sequences.
## Processing: NTC_13 
## Sample 1 - 3348 reads in 1724 unique sequences.
## Sample 1 - 3348 reads in 1565 unique sequences.
## Processing: NTC_14 
## Sample 1 - 1466 reads in 631 unique sequences.
## Sample 1 - 1466 reads in 605 unique sequences.
## Processing: NTC_15 
## Sample 1 - 355 reads in 242 unique sequences.
## Sample 1 - 355 reads in 235 unique sequences.
## Processing: NTC_16 
## Sample 1 - 934 reads in 517 unique sequences.
## Sample 1 - 934 reads in 475 unique sequences.
## Processing: NTC_17 
## Sample 1 - 1450 reads in 951 unique sequences.
## Sample 1 - 1450 reads in 852 unique sequences.
## Processing: NTC_18 
## Sample 1 - 294 reads in 202 unique sequences.
## Sample 1 - 294 reads in 199 unique sequences.
## Processing: NTC_19 
## Sample 1 - 782 reads in 740 unique sequences.
## Sample 1 - 782 reads in 734 unique sequences.
## Processing: NTC_2 
## Sample 1 - 4900 reads in 1632 unique sequences.
## Sample 1 - 4900 reads in 1215 unique sequences.
## Processing: NTC_20 
## Sample 1 - 874 reads in 575 unique sequences.
## Sample 1 - 874 reads in 580 unique sequences.
## Processing: NTC_3 
## Sample 1 - 2041 reads in 1064 unique sequences.
## Sample 1 - 2041 reads in 904 unique sequences.
## Processing: NTC_4 
## Sample 1 - 305 reads in 160 unique sequences.
## Sample 1 - 305 reads in 160 unique sequences.
## Processing: NTC_5 
## Sample 1 - 933 reads in 328 unique sequences.
## Sample 1 - 933 reads in 231 unique sequences.
## Processing: NTC_6 
## Sample 1 - 1023 reads in 366 unique sequences.
## Sample 1 - 1023 reads in 257 unique sequences.
## Processing: NTC_7 
## Sample 1 - 1262 reads in 574 unique sequences.
## Sample 1 - 1262 reads in 571 unique sequences.
## Processing: NTC_8 
## Sample 1 - 1286 reads in 702 unique sequences.
## Sample 1 - 1286 reads in 593 unique sequences.
## Processing: NTC_9 
## Sample 1 - 196 reads in 145 unique sequences.
## Sample 1 - 196 reads in 136 unique sequences.
## Processing: OM18_BC 
## Sample 1 - 4691 reads in 1345 unique sequences.
## Sample 1 - 4691 reads in 1132 unique sequences.
## Processing: OM18_BD 
## Sample 1 - 2620 reads in 1606 unique sequences.
## Sample 1 - 2620 reads in 1587 unique sequences.
## Processing: OM18_BJ 
## Sample 1 - 1483 reads in 1221 unique sequences.
## Sample 1 - 1483 reads in 1200 unique sequences.
## Processing: WAB103_10_10D 
## Sample 1 - 22517 reads in 7639 unique sequences.
## Sample 1 - 22517 reads in 5308 unique sequences.
## Processing: WAB103_22_10C 
## Sample 1 - 29843 reads in 10230 unique sequences.
## Sample 1 - 29843 reads in 9504 unique sequences.
## Processing: WAB103_45_10J 
## Sample 1 - 29545 reads in 12048 unique sequences.
## Sample 1 - 29545 reads in 10343 unique sequences.
## Processing: WAB104_10_5J 
## Sample 1 - 35972 reads in 12607 unique sequences.
## Sample 1 - 35972 reads in 10735 unique sequences.
## Processing: WAB104_22_5D 
## Sample 1 - 26450 reads in 7145 unique sequences.
## Sample 1 - 26450 reads in 5931 unique sequences.
## Processing: WAB104_45_5C 
## Sample 1 - 41957 reads in 15201 unique sequences.
## Sample 1 - 41957 reads in 12948 unique sequences.
## Processing: WAB105_10_6C 
## Sample 1 - 30423 reads in 6375 unique sequences.
## Sample 1 - 30423 reads in 5315 unique sequences.
## Processing: WAB105_22_6J 
## Sample 1 - 43572 reads in 9187 unique sequences.
## Sample 1 - 43572 reads in 7875 unique sequences.
## Processing: WAB105_45_6D 
## Sample 1 - 36153 reads in 7575 unique sequences.
## Sample 1 - 36153 reads in 6784 unique sequences.
## Processing: WAB188_10_4D 
## Sample 1 - 29932 reads in 6019 unique sequences.
## Sample 1 - 29932 reads in 4747 unique sequences.
## Processing: WAB188_22_4C 
## Sample 1 - 13245 reads in 4398 unique sequences.
## Sample 1 - 13245 reads in 3855 unique sequences.
## Processing: WAB188_45_4J 
## Sample 1 - 37216 reads in 12296 unique sequences.
## Sample 1 - 37216 reads in 10561 unique sequences.
## Processing: WAB55_10_8J 
## Sample 1 - 20226 reads in 6186 unique sequences.
## Sample 1 - 20226 reads in 5134 unique sequences.
## Processing: WAB55_22_8D 
## Sample 1 - 44540 reads in 13655 unique sequences.
## Sample 1 - 44540 reads in 11528 unique sequences.
## Processing: WAB55_45_8C 
## Sample 1 - 29464 reads in 9605 unique sequences.
## Sample 1 - 29464 reads in 7856 unique sequences.
## Processing: WAB71_10_3C 
## Sample 1 - 22978 reads in 6253 unique sequences.
## Sample 1 - 22978 reads in 6052 unique sequences.
## Processing: WAB71_22_3J 
## Sample 1 - 30053 reads in 6642 unique sequences.
## Sample 1 - 30053 reads in 5972 unique sequences.
## Processing: WAB71_45_3D 
## Sample 1 - 41879 reads in 10131 unique sequences.
## Sample 1 - 41879 reads in 9224 unique sequences.
rm(derepF); rm(derepR)

6.3 Construct sequence table

seqtab <- makeSequenceTable(mergers)
write_rds(seqtab, path = format(Sys.Date(), "data_output/seqtab_OM18_processed_%Y%m%d.rds"), compress = "gz")

7 Remove Chimeras

Although dada2 has searched for insertion/deletion errors and substitutions, there may still be chimeric sequences in our dataset (sequences that are derived from forward and reverse sequences from two different organisms becoming fused together during PCR and/or sequencing). To identify chimeras, we will search for rare sequence variants that can be reconstructed by combining left-hand and right-hand segments from two more abundant “parent” sequences.

# Remove chimeras
seqtab_nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE)
# Print percentage of our sequences that were not chimeric.
100*sum(seqtab_nochim)/sum(seqtab)
## [1] 98.35191
write_rds(seqtab_nochim, path = format(Sys.Date(), "data_output/seqtab_nochim_OM18_processed_%Y%m%d.rds"), compress = "gz")

8 Track reads through pipeline

getN <- function(x) sum(getUniques(x)) # function to grab sequence counts from output objects
# tracking reads by counts
filt_out_track <- filt_out %>%
  data.frame() %>%
  mutate(Sample = rownames(.) %>% str_remove(pattern = "_S([[:digit:]]+)_L001_R1_001.fastq.gz")) %>%
  rename(input = reads.in, filtered = reads.out)
rownames(filt_out_track) <- filt_out_track$Sample

ddF_track <- data.frame(denoisedF = sapply(ddF[sample_names_F_filt], getN)) %>%
  mutate(Sample = row.names(.))
ddR_track <- data.frame(denoisedR = sapply(ddR[sample_names_F_filt], getN)) %>%
  mutate(Sample = row.names(.))
merge_track <- data.frame(merged = sapply(mergers, getN)) %>%
  mutate(Sample = row.names(.))
chim_track <- data.frame(nonchim = rowSums(seqtab_nochim)) %>%
  mutate(Sample = row.names(.))

track <- left_join(filt_out_track, ddF_track, by = "Sample") %>%
  # left_join(ddF_track, by = "Sample") %>%
  left_join(ddR_track, by = "Sample") %>%
  left_join(merge_track, by = "Sample") %>%
  left_join(chim_track, by = "Sample") %>%
  replace(., is.na(.), 0) %>%
  select(Sample, everything())
row.names(track) <- track$Sample
head(track)
# tracking reads by percentage
track_pct <- track %>% 
  data.frame() %>%
  mutate(Sample = rownames(.),
         filtered_pct = ifelse(filtered == 0, 0, 100 * (filtered/input)),
          # cut_pct = ifelse(cut == 0, 0, 100 * (cut/filtered)),
         denoisedF_pct = ifelse(denoisedF == 0, 0, 100 * (denoisedF/filtered)),
         denoisedR_pct = ifelse(denoisedR == 0, 0, 100 * (denoisedR/filtered)),
         merged_pct = ifelse(merged == 0, 0, 100 * merged/((denoisedF + denoisedR)/2)),
         nonchim_pct = ifelse(nonchim == 0, 0, 100 * (nonchim/merged)),
         total_pct = ifelse(nonchim == 0, 0, 100 * nonchim/input)) %>%
  select(Sample, ends_with("_pct"))

# summary stats of tracked reads averaged across samples
track_pct_avg <- track_pct %>% summarize_at(vars(ends_with("_pct")), 
                                            list(avg = mean))
head(track_pct_avg)
track_pct_med <- track_pct %>% summarize_at(vars(ends_with("_pct")), 
                                            list(avg = stats::median))
head(track_pct_avg)
head(track_pct_med)
# Plotting each sample's reads through the pipeline
track_plot <- track %>% 
  pivot_longer(-Sample, names_to = "Step", values_to = "Reads") %>% 
  mutate(Step = factor(Step, 
                       levels = c("input", "filtered", "denoisedF", "denoisedR", "merged", "nonchim"))) %>%
  ggplot(aes(x = Step, y = Reads)) +
  geom_line(aes(group = Sample), alpha = 0.2) +
  geom_point(alpha = 0.5, position = position_jitter(width = 0)) + 
  stat_summary(fun = median, geom = "line", group = 1, color = "steelblue", size = 1, alpha = 0.5) +
  stat_summary(fun = median, geom = "point", group = 1, color = "steelblue", size = 2, alpha = 0.5) +
  stat_summary(fun.data = median_hilow, fun.args = list(conf.int = 0.5), 
               geom = "ribbon", group = 1, fill = "steelblue", alpha = 0.2) +
  geom_label(data = t(track_pct_avg[1:5]) %>% data.frame() %>% 
               rename(Percent = 1) %>%
               mutate(Step = c("filtered", "denoisedF", "denoisedR", "merged", "nonchim"),
                      Percent = paste(round(Percent, 2), "%")),
             aes(label = Percent), y = 1.1 * max(track[,2])) +
  geom_label(data = track_pct_avg[6] %>% data.frame() %>%
               rename(total = 1),
             aes(label = paste("Total\nRemaining:\n", round(track_pct_avg[1,6], 2), "%")), 
             y = mean(track[,6]), x = 6.6) +
  expand_limits(y = 1.1 * max(track[,2]), x = 7.2) +
  theme_classic()

track_plot